def read_particles_info(filename):
    particles = []
    with open(filename, 'r') as file:
        for line in file:
            data = line.split()
            N = int(data[0])
            Q = float(data[1])
            x, y, z = map(float, data[2:])
            particles.append((N, Q, x, y, z))
    return particles


def read_mesh_info(filename):
    nodes = []
    elements = []
    with open(filename, 'r') as file:
        num_nodes = int(file.readline())
        for _ in range(num_nodes):
            data = file.readline().split()
            x, y, z = map(float, data)
            nodes.append((x, y, z))

        num_elements = int(file.readline())
        for _ in range(num_elements):
            data = file.readline().split()
            n1, n2, n3, n4 = map(int, data)
            elements.append((n1, n2, n3, n4))
    return nodes, elements


def compute_average_charge_density(particles, nodes, elements):
    densities = [0] * len(elements)
    for particle in particles:
        N, Q, x, y, z = particle
        for i, element in enumerate(elements):
            n1, n2, n3, n4 = element
            x1, y1, z1 = nodes[n1 - 1]
            x2, y2, z2 = nodes[n2 - 1]
            x3, y3, z3 = nodes[n3 - 1]
            x4, y4, z4 = nodes[n4 - 1]

            volume = compute_tetrahedron_volume(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4)
            density = Q / volume
            densities[i] += density

    average_densities = [density / len(particles) for density in densities]
    return average_densities


def compute_tetrahedron_volume(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4):
    volume = (1 / 6) * abs((x2 - x1) * ((y3 - y1) * (z4 - z1) - (z3 - z1) * (y4 - y1))
                           + (y2 - y1) * ((z3 - z1) * (x4 - x1) - (x3 - x1) * (z4 - z1))
                           + (z2 - z1) * ((x3 - x1) * (y4 - y1) - (y3 - y1) * (x4 - x1)))
    return volume


def find_element_by_xyz(x, y, z, nodes, elements):
    for i, element in enumerate(elements):
        n1, n2, n3, n4 = element
        x1, y1, z1 = nodes[n1 - 1]
        x2, y2, z2 = nodes[n2 - 1]
        x3, y3, z3 = nodes[n3 - 1]
        x4, y4, z4 = nodes[n4 - 1]

        if is_point_inside_tetrahedron(x, y, z, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4):
            return i

    return -1


def is_point_inside_tetrahedron(x, y, z, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4):
    v1 = (x1 - x, y1 - y, z1 - z)
    v2 = (x2 - x, y2 - y, z2 - z)
    v3 = (x3 - x, y3 - y, z3 - z)
    v4 = (x4 - x, y4 - y, z4 - z)

    dot_product1 = dot_product(v1, cross_product(v2, v3))
    dot_product2 = dot_product(v2, cross_product(v3, v4))
    dot_product3 = dot_product(v3, cross_product(v4, v1))
    dot_product4 = dot_product(v4, cross_product(v1, v2))

    return dot_product1 >= 0 and dot_product2 >= 0 and dot_product3 >= 0 and dot_product4 >= 0


def cross_product(v1, v2):
    x = v1[1] * v2[2] - v1[2] * v2[1]
    y = v1[2] * v2[0] - v1[0] * v2[2]
    z = v1[0] * v2[1] - v1[1] * v2[0]
    return x, y, z


def dot_product(v1, v2):
    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]


# 读取粒子信息文件
particles = read_particles_info('particles.txt')

# 读取网格文件
nodes, elements = read_mesh_info('mesh.txt')

# 统计每个有限元的平均电荷密度
average_densities = compute_average_charge_density(particles, nodes, elements)
print('Average Charge Densities:')
for i, density in enumerate(average_densities):
    print(f'Element {i + 1}: {density}')

# 输入外部的xyz，判断所在有限元并输出对应的平均电荷密度值
x = float(input('Enter x: '))
y = float(input('Enter y: '))
z = float(input('Enter z: '))

element_index = find_element_by_xyz(x, y, z, nodes, elements)
if element_index != -1:
    average_density = average_densities[element_index]
    print(f'Average Charge Density at ({x}, {y}, {z}): {average_density}')
else:
    print('The point is not inside any element.')

# 请将粒子信息文件命名为particles.txt，将网格文件命名为mesh.txt，
# 并确保它们与脚本文件在同一目录中。运行脚本后，将按照要求进行操作并得到相应的输出。
